Modelling liver cancer microenvironment using a novel 3D culture system

The tumor microenvironment and its contribution to tumorigenesis has been a focal highlight in recent years. A two-way communication between the tumor and the surrounding microenvironment sustains and contributes to the growth and metastasis of tumors. Progression and metastasis of hepatocellular carcinoma (HCC) have been reported to be exceedingly influenced by diverse microenvironmental cues. In this study, we present a 3D-culture model of liver cancer to better mimic in vivo tumor settings. By creating novel 3D co-culture model that combines free-floating and scaffold-based 3D-culture techniques of liver cancer cells and fibroblasts, we aimed to establish a simple albeit reproducible ex vivo cancer microenvironment model that captures tumor-stroma interactions. The model presented herein exhibited unique gene expression and protein expression profiles when compared to 2D and 3D mono-cultures of liver cancer cells. Our results showed that in vivo like conditions cannot be mimicked by simply growing cancer cells as spheroids, but by co-culturing them with 3D fibroblast with which they were able to crosstalk. This was evident by the upregulation of several pathways involved in HCC, and the increase in secreted factors by co-cultured cancer cells, many of which are also involved in tumor-stroma interactions. Compared to the conventional 2D culture, the proposed model exhibits an increase in the expression of genes associated with development, progression, and poor prognosis of HCC. Our results correlated with an aggressive outcome that better mirrors in vivo HCC, and therefore, a more reliable platform for molecular understanding of HCC.

A tumor has an increasing demand for oxygen and nutrients to support its progression. When the demand for oxygen remains unmet, low oxygen hypoxic conditions ensue 12 . To survive, tumor cells activate the hypoxiainducible factor 1 (HIF1) 12 , which in turn activates the transcription of a group of genes through binding to their hypoxia-response elements to promote the survival of tumor cells 13 . HIF-1 targeted genes significantly contribute to tumor angiogenesis, metastasis, adhesion, metabolism, and pH regulation 13 . Moreover, many studies have highlighted the role of hypoxia in recruiting stromal components to the TME 14 , ECM composition, and metastatic remodeling 15 .
Hepatocellular carcinoma (HCC), is the fifth most common cancer and is the fourth cause of cancer-related death worldwide 16 . HCC has a very poor prognosis with only five-year survival rate 17 . HCC progression is influenced by the liver microenvironment such as altered stromal cells 18 . These cells deposit ECM proteins causing fibrosis that then progresses to cirrhosis 19 ; suggesting a crucial role of ECM build-up in HCC progression 20 . Hypoxia represents a driving force for HCC progression, and is associated with poor prognosis 21 . HIF-mediated gene expression contributes to different aspects of HCC metastasis, such as epithelial mesenchymal transition (EMT) 22 , invasion of the ECM, and metastasis 23 . Yet the molecular mechanisms governing stromal and tumor cell interactions within the TME of HCC under hypoxic conditions 14 remains unclear.
To reflect the complexity and dynamic nature of tumor cell biology, a physiologically relevant model is needed. Especially, when it comes to drug discovery and identifying effective therapeutic targets. To simulate in vivo environment, in vitro two-dimensional (2D) cell culture is typically assembled by growing cells on a plastic substrate as an adherent monolayer. However, distortion of spatial arrangement of cells in 2D culture changes cell-cell and cell-matrix interactions 24 , and most importantly, alters the response of cells to certain drugs and treatments 25 . That is why, 3D cultured cells better recapitulates in vivo architecture of tumors and exhibits gene expression closer to that of in vivo tumors 24 . One very common 3D culture cell model is the spheroid; a micro cell cluster sphere 26 . The nature of this three-dimensional multicellular model is what makes it an attractive tool to simulate solid tumors in vitro as it is composed of three regions, a highly proliferative outer region, a middle quiescent region, and a hypoxic core region 27 . Such compartmentalization creates diffusional gradients of oxygen, nutrients, and tested drugs among all three regions of the spheroid, which is also characteristic of solid tumors 28 .
We aimed at modeling the basic TME of liver cancer by mimicking certain aspects of in vivo tumors, such as three-dimensionality of tissue, hypoxia, and heterogeneity of tumors. Five groups were designed to reflect each element, group 1 is a control group for comparison purposes which consists of 2D mono-cultures of liver cancer cells. Group 2 is also 2D mono-cultures of liver cancer cells but chemically induced for hypoxia. Group 3 on the other hand is like group 2 but additionally includes conditioned media from 2D fibroblasts to reflect a one-way co-culture system. Group 4 and 5 are 3D cultures of liver cancer cells that exhibit hypoxia physiologically due to culturing conditions. However, group 5, which is our proposed model, includes a 3D culture of fibroblast in a separate insert, reflecting a two-way co-culture system. Our findings reported herein demonstrate that our proposed model of group 5 reflects many aspects of in vivo settings and signaling pathways, promoting it as a potential platform for further studies of drug efficacy in vitro and understanding the communication between cancer and the stroma in liver cancer.

Materials and methods
Cell lines and cell culture. HepG2 (human Hepatoma cell line) and SV-80 (human fibroblast) were purchased from CLS (CLS GmbH, Germany). HepG2 and SV-80 cells were maintained in high-glucose DMEM medium supplemented with 10% FBS and 1% antibiotic/ antimycotic cocktail (HyClone, UK) at 37 °C in 5% CO 2 humidified incubator. Cells were sub-cultured every 2-4 days using trypsin 0.25%. For 2D mono-cultures, HepG2 and SV-80 cells were cultured at a density of 1 × 10 6 in conventional 2D culturing flasks in DMEM media. Flasks were incubated at 37 °C in 5% CO 2 humidified incubator. For generation of HepG2 spheroids, HepG2 cells were trypsinized and resuspended as single cell suspension before being seeded at 1 × 10 6 cell density in Corning® Ultra-Low attachment cell culture flasks coated with poly-HEMA (Corning, USA), generating spheroids of different sizes ( Fig. S1a-b). Plates were incubated at 37 °C in 5% CO 2 humidified incubator and spheroid formation was monitored over a period of 5 days using inverted microscopy. For generation of 3D culture of fibroblasts, Alvetex strata inserts (Reinnervate, UK) were used. SV-80 cells were trypsinized and resuspended as single cell suspension for cell counting. Alvetex inserts were prepared prior to cell seeding by three washes (1st: 70% ethanol, 2nd: growth media, 3rd: growth media). Inserts were placed in 6-well plate and SV-80 cells were seeded in the Alvetex inserts in a density of 1 × 10 6 .

Morphology assessment.
To study the morphological differences between 2D and 3D cultures, HepG2 and SV-80 cells were seeded in in 8-well chambers for 2D culture; or in ULA plates and Alvetex inserts, respectively, for 3D culture. 2D cultures of HepG2 and SV-80 were seeded at a density of 2 × 10 4 and 1 × 10 4 cells/ well in an 8-well chamber, respectively. They were fixed with 100% methanol and stained with crystal violet. 3D culture of HepG2 was seeded at a density of 1 × 10 6 and monitored over a period of 5 days. Spheroids were harvested, fixed, and stained with crystal violet and immobilized on agarose pads for imaging. 3D culture of SV-80 was seeded at a density of 0.5 × 10 6 , then fixed and stained with neutral red. Inserts were unclipped and scaffolds were placed on glass slides for imaging using IX53 inverted microscope (Olympus, Japan).
Hypoxia-mimicking conditions. To mimic the hypoxic microenvironment, 2D cultures of HepG2 cells were cultured as described above and treated with 200 and 300 μM of Cobalt (II) Chloride hexahydrate (CoCl 2 ) for 6 h prior to harvesting. The doses and incubation time are based on literature demonstrating HIF1-α maximum induction at 4-6 h of CoCl 2 treatment, over a range of doses. Under normoxia, HIF1-α protein is degraded, Immunofluorescence. HepG2 cells were seeded at a density of 1 × 10 4 cells/ well in an 8-well chamber and allowed to reach confluency. Cells were incubated with or without CoCl 2 for 6 h, after which they were washed with ice-cold PBS and then fixed with cold absolute methanol for 10 min at − 20 °C. Cells were washed and incubated with 1% BSA blocking solution for 30 min at room temperature. Blocking solution was discarded and cells were washed and incubated with antibody against HIF1α (ab1) at 1:50 dilution at 4 °C overnight. Cells were washed and incubated with secondary antibody prior to counter-staining with DAPI and imaging.
Co-culture systems. For the 2D co-culture, a one-way communication system was followed. Briefly, HepG2 and SV-80 cells were seeded at a density of 1 × 10 6 in conventional 2D culturing flasks in DMEM media. Flasks were incubated at 37 °C in 5% CO 2 humidified incubator for about 24 h. Conditioned media of SV-80 cell line was collected, centrifuged to collect any cellular debris, and applied to HepG2 cells, which were incubated with the conditioned media for 48 h.
For the 3D co-culture, a two-way communication system was followed. Briefly, HepG2 and SV-80 3D cultures were prepared separately. Prior to co-culture, 6-well plates were coated with 1.5% agarose, and allowed to set and cool before transferring HepG2 spheres to the bottom of the coated plates, and the inserts containing SV-80 3D culture were placed on top. HepG2 and SV-80 3D cultures were incubated for 48 h.

RNA-seq libraries construction and sequencing.
For RNA extraction, all groups were prepared as described previously, then collected, washed with 1X PBS, and resuspended in RNAlater stabilization solution before storing at − 80 °C. Total RNA was isolated from three biological replicates of all groups using RNeasy Mini Kit (Qiagen) following manufacturer's instructions. Concentration and purity of total RNA was assessed using NanoDrop2000. Quality control of RNA samples was performed with Agilent 2100 Bioanalyzer RNA 6000 Nano Kit, concentrated samples were diluted with RNase-free water prior to bioanalyzer run. The RNAseq libraries were prepared by DNA Sequencing Center in Brigham Young University. Briefly, KAPA Stranded mRNA-Seq Kit (Kapa Biosystems, USA) was used for capturing poly(A) RNA, converting it to cDNA, A-tailing, and Adapter ligation. Fragments carrying appropriate adapter sequences were amplified to yield mRNA-Seq libraries. KAPA Library Quantification Kit (Kapa Biosystems, USA) was used for libraries quantification prior to Illumina sequencing using high-throughput Illumina HiSeq sequencing system (Illumina, USA).

Alignment and analysis of illumina reads.
Briefly, reads obtained from Illumina were aligned to Homo sapiens GRCh38.p2 reference genome using tophat2 v2.1.0. Following alignment and annotation, read counts were generated using HTseq count. Read counts were used to generate principle component analysis (PCA) plot and hierarchical cluster heatmaps, using the web-based tool ClustVis (https:// biit. cs. ut. ee/ clust vis_ large/), for clustering of multivariate data. Triplicates of each group are collapsed by taking the mean, and rows were scaled using vector scaling method.
Quantitative polymerase chain reaction. For the purpose of validating genes of interest, quantitative polymerase chain reaction (qPCR) was carried out on samples from groups 1 and 5. Total RNA was converted to cDNA using GoScript™ Reverse Transcription System, according to manufacturer's instruction (Promega, USA). GoTaq® qPCR Master Mix (Promega) was used to perform qPCR on QuantStudio 5 Real-Time PCR System (Applied Biosystems) using gene-specific primers purchased from Macrogen (Macrogen Inc.). 18S rRNA was used for data normalization due to its unchanged expression in groups 1 and 5. Comparative C T method ( www.nature.com/scientificreports/ was used to determine fold change in expression of target genes between group 5 and group 1, according to the following equation: 2 -ΔΔC T = [(C T gene of interest-C T internal control) group 5-(C T gene of interest-C T internal control) group 1]. The cycling parameters recommended in GoTaq® qPCR Master Mix manual were used. Briefly, Hot-Start Activation was carried at 95 °C for 2 min, 40 cycles of denaturation at 95 °C for 15 s, followed by annealing/extension at 60 °C for 60 s, and then dissociation at 60-95 °C.
Gene set and gene ontology enrichment analyses. Gene set enrichment analysis (GSEA) was carried on sets of upregulated and downregulated genes using canonical pathways ontology on eXploring Genomic Relations (XGR) web version (http:// galah ad. well. ox. ac. uk: 3030/). Enriched terms were tested for significance using the Hypergeometric test, and only terms with FDR < 0.05 were considered. Gene ontology (GO) enrichment analysis was carried out using Biological Networks Gene Ontology (BiNGO) App 29 in Cytoscape, an open source software platform for network data integration, analysis, and visualization 30 . Based on Hypergeometric significance test, corrected multiple testing using Benjamini and Hochberg FDR < 0.05, and using biological process ontology, relevant enriched terms were identified by BiNGO. InteractiVenn (http:// www. inter activ enn. net/) was used to generate a Venn diagram of relevant GO terms by uploading sets of their associated genes 31  Functional analysis. Functional analyses were performed on differential expressed genes with FDR < 0.05 from the different group comparisons (G2 vs G1, G3 vs G1, G4 vs G1 and G5 vs G1). Gene onthology (GO) and pathway clustering analysis were performed using the ClueGO, Cytoscape plug-in; pathway analyses were based on the Kyoto encyclopedia of genes and genomes (KEGG) [32][33][34] and WikiPathways annotations built-in ClueGO 35 . CluePedia Cytoscape plug-in, was used on ClueGO outcome to show the molecules tighten to the highlighted nodes/networks 36 . ClueGO/CluePedia analyses followed default parameters (or otherwise in-text specified). Term convergence among GO, KEGG and Wikipathways terms assures the validity of the model created to each of the group comparisons.
Pathway analysis. Pathway visualization and analysis were performed in PathVisio 3.3.0. These analyses were based on the Wikipathways human collection. Primary analyses were carried out in order to visualize how the most significant molecules included in ClueGO/CluePedia analysis fit on the most significant pathways identified through this functional clustering. To import this data into Pathvisio, as an input the following data were included per molecule in a table format: EnsEMBL ID, FDR and logFC. After pathway mapping and identification of these molecules within the pathway maps of interest, all the molecules with an FDR below 0.05 were included into Pathvisio in order to fill the gaps within the focused pathways and proceed to identification of trends in regulation. LogFC values from 1 to − 1 are represented and visualized in a color gradient from red to green respectively. Kaplan-Meier survival analysis were performed on GDC TCGA Liver Cancer (LIHC) primary tumor samples (n = 179) for genes within our main pathway of interest (Edited on Pathvisio) through UCSC Xena browser 37 . Kaplan-Meier curves are compared using the log-rank test. Here, the UCSC Xena browser reported p-value (χ 2 distribution) is shown.
Growth factor antibody array. Secreted growth factors, angiogenesis factors, and cytokines in conditioned media of 3D mono-or co-cultures were analyzed using Human Growth Factor Antibody Array kit (ab134002, abcam), human Angiogenesis Antibody Array kits (ab169808/ab134000, abcam). Briefly, membranes were incubated with the blocking solution for 30 min at room temperature. After which membranes were washed and incubated with 2 mL of conditioned media collected from indicated groups overnight at 4 °C. Membranes were washed and incubated with Biotin-Conjugated Anti-Cytokines overnight at 4 °C. The membranes were the incubated with 2 mL of HRP-Conjugated Streptavidin for 2 h at room temperature. Membranes were washed prior to incubation with detection buffer for 2 min and imaged using ChemiDoc imaging system (Biorad). Mean pixel intensity were determined using ImageJ.
Statistical analysis. Statistical analysis was performed using Microsoft Excel. Two-tailed, paired t-test was used to analyze the data, where statistical significance was assumed at *p < 0.01 **p < 0.001. Data is represented as mean ± SD.

Results
Comparing experimental setup of 2D, mono-3D, and co-3D cultures. Fibroblasts grown as a 2D culture exhibit an elongated morphology. However, growing fibroblasts as 3D culture in porous scaffolds alters their morphology to be more rounded. HepG2 cells grown in 3D were monitored over a period of 5 days to assess formation of tight spheroids with a smooth surface (Fig. 1a,b) www.nature.com/scientificreports/ increasing concentrations of CoCl 2 (100-400 μM) to assess cellular viability under hypoxia-mimicking conditions (Fig. 1c). Treatment of HepG2 with CoCl 2 did not affect cellular viability significantly at doses of 100 and 200 μM of CoCl 2 . However, a highly significant difference (p < 0.001) was noted at a dose of 400 μM CoCl 2 .
To confirm the induction of hypoxia in CoCl 2 -treated HepG2 2D cultures, protein expression of HIF1α was assessed using western blot. Treatment with 200 μM of CoCl 2 for 6 h did not induce HIF1α expression. However, by increasing the dose to 300 μM, HIF1α expression was detected in HepG2 cells (Fig. 1d). Immunofluorescence was then used to detect the cellular localization of HIF1α, showing that CoCl 2 also affects HIF1α translocation to the nucleus (Fig. 1e), where it can bind to hypoxia-response elements (HREs).
Differential gene expression analysis. To better assess in vivo mimicking capabilities of our 5 groups, RNA from the denoted samples ( Fig. 1a) was extracted and then sequenced via the Illumina platform. Significant differentially expressed genes (DEGs) were determined at a cutoff of Log2FC of 1 and FDR < 0.05. A hierarchical clustering Heatmap was used to explore changes in global gene expression across our 5 models (Fig. 2a).
Interestingly, group 5 clustered separately from the rest, revealing a fundamental and distinctive change in gene expression. Principal component analysis (PCA) based on the global gene expression showed a similar outcome, where group 5 segregated at the opposite extreme of all other groups, with the conventional HepG2 culture (group 1) being the most distant (Fig. 2b), suggesting fundamental differences between 2D, monocellular 3D, and multicellular 3D models. Significant differentially expressed genes (DEG) were determined relative to group 1 using RNA-seq 2G. Group2 resulted in significant upregulation of 243 genes and downregulation of 131 genes. Group3 increased the number of significantly DEGs to 474 upregulated and 145 downregulated genes. When comparing gene  (Fig. 2c). In Addition, significant DEGs signatures, unique to each group, were identified (Fig. 2d). Group 5 exhibited the highest number in unique significant DEGs among all other groups, at 957 upregulated genes and 796 downregulated genes. A full List of common and unique significantly DEGs of each group is available in supporting files 1 and 2. Canonical pathways associated with significant DEGs in different culturing conditions were analyzed using gene-based enrichment analysis by XGR. As expected, culturing HepG2 cells under hypoxia-mimicking www.nature.com/scientificreports/ conditions (group 2 and 3) upregulated genes involved in hypoxia-inducible factor-1 alpha (HIF1-α) and hypoxia-inducible factor-2 alpha (HIF2-α) pathways, and networks downstream of these pathways. This was even observed in 3D spheroid cultures (group 4 and 5) despite not being treated with a hypoxia inducing agent, indicating the formation of hypoxic core in 3D culture spheroids. Culturing HepG2 cells with only fibroblasts conditioned media (group 3) or alternatively with 3D culture of fibroblasts (group 5) upregulated genes involved in integrin family cell surface interactions, interleukin-6 (IL6) mediated signaling events (Table S1).

In-depth pathway functional analysis reveals pathways associated with HCC progression.
To further understand the role of significant DEGs in each group, we applied clustering analysis using ClueGO/ CluePedia as described previously 38 . By combining Gene Ontology (GO) terms, KEGG and Wiki pathways, ClueGo/CluePedia create a better interpretation of the pathways associated with the list of input genes 35 . Hypoxia mimicking conditions in group 2 induced cellular responses to hypoxia and HIF-1 signaling pathway, as well as other signaling pathways reported to promote HCC including NRF2, FOXO, and p53 pathways 39,40 (Fig. 3a,b). Similarly, inducing hypoxia in group 3 resulted in inducing HIF-1 signaling pathways. However, with the addition of fibroblast conditioned media, group 3 DEGs-associated processes were enriched in 3 out of 5 pathways previously reported in KEGG analysis of HCC patients' tissues including (i) complement and coagulation cascades, (ii) focal adhesion, and (iii) ECM-receptor interaction 41 (Fig. 3a). DEGs-associated processes of HepG2 3D culture alone (group 4) resembled some of group 2 such as hypoxia and NRF2 signaling pathways, but also some of group 3 such as G3 such as complement and coagulation cascades (Fig. 3a). Additionally, DEGsassociated processes of group 4 included steroid hormone biosynthesis process and estrogen signaling pathway, which have been linked to HCC progression 42 (Fig. 3a). Group 5 shared some pathways with the other groups such as hypoxia, focal adhesion, glycolysis/gluconeogenesis, and estrogen signaling pathway. In addition, DEGsassociated processes of group 5 were significantly enriched in HCC-promoting pathways including oncostatin M signaling pathway, insulin signaling pathway and aryl hydrocarbon receptor pathways 43,44 (Fig. 3b).
After identifying the main pathways involved in group 5 through ClueGO/CluePedia clustering, Insulin signaling pathway was further analyzed and visualized in PathVisio due to its documented relevance in HCC 45 and to gain insights of significant physiological changes occurring within a specific pathway. Genes below FDR 0.05 (5871 hits) were imported into PathVisio 46 to identify trends in regulation within this particular pathway and other chained events embedded within associated pathways. Inconsistencies within these maps were excluded from the final pathway (Fig. S4). The Insulin signaling pathway in group 5 was found to mainly lead to activation of genes involved MAPK signaling and this trend converges with hypoxia signaling input, both leading to the production of VEGFB, VEGFA, PGF, growth factors known to be involved in angiogenesis and tumor invasion 47 . Tumor invasion can also be promoted via Ras and its downstream pathways of MAPK and RALB, which were all shown to be upregulated in group 5. Some of the genes involved in the pathway shown previously have been validated with qPCR (Fig. 4a).

Co-culture of fibroblasts and HepG2 3D cultures enriches the secretome in pro-tumor factors.
To identify signaling factors playing a role in the crosstalk between cancer cells and fibroblasts, the secretome of group 4 and 5 was analyzed for factors involved in angiogenesis, invasion, and metastasis, in addition to other factors in the insulin signaling pathway, secretome of 3D fibroblasts was included as a control (Fig S2, Fig. 4b-f). Most factors were dramatically increased in the secretome of group 5 (from both fibroblasts and HepG2 cells) when compared to 3D mono-cultures. Co-culturing HepG2 spheroids with 3D culture of fibroblasts increased the levels of insulin signaling pathway factors like IGF-II, IGFBP1, IGFBP2 (Fig. 4b). In addition, levels of different angiogenesis and cytokines that are involved in the crosstalk in the TME were highly increased in the setting of group 5 (Fig. 4c,f).
As the secretome collected from group 5 is shared among fibroblasts and HepG2 cells, we wanted to confirm the source of the secreted factors. To that end we compared gene expression of 3D HepG2 or 3D fibroblasts before and after co-culture (group 5 over group 4) (Fig. 4g). Differential gene expression analysis revealed that M-CSF, IGFBP1, IGFBP2, TGF-β, and PECAM-1 were upregulated in HepG2 cells after co-culture, while they were not differentially expressed in fibroblasts. On the other hand, TGF-α, and IL-1β were only differentially expressed in fibroblasts after co-culture. VEGF and IGF-II were differentially expressed in both HepG2 cells and fibroblasts, being more upregulated in HepG2 cells (Fig. 4g). Despite that VEGF is upregulated more than 2.5 folds on the RNA level in HepG2 cells after co-culture (Fig. 4g), there was no noticeable difference in the secreted VEGF from 3D HepG2 cells before and after co-culture (Fig. 4f). Nonetheless, VEGFR2 was more expressed in group 5 in comparison to group 1 and group 4; while VEGFR1 was higher in the 3D cultures, but no noticeable difference between the two (Fig. 4h). TNFR1 was also dramatically higher in group 5 when compared to both group 1 and group 4 (Fig. 4h). Kaplan-Meier curves of genes within our main pathway of interest were analyzed. IGF1R and EGLN3 showed a significant correlation between expression and prognosis (Fig. 4i,j). Both genes were shown to be involved in insulin signaling pathway (Fig. S4).

Pathway complementation analysis highlights regulatory miRNAs and prognosis markers.
To better understand the changes elicited by the 3D co-culture conditions on HepG2 cells, we investigated the effect of miRNA regulation given their role in affecting de novo or modulating established gene expression in tumors 48,49 . Based on miRNA-gene interaction analysis of significant DEGs in group 5, miR-335 was the top regulatory miRNA with the highest number of connections to other genes in the network (Fig. S5a). Other miRNAs were predicted, top 10 are shown in Fig. S5b  www.nature.com/scientificreports/  [32][33][34] . All groups are compared to G1 as the control condition. Node size of pathway terms resemble the number of associated genes to it. The stronger the node color the more significant a cluster is. (b) ORA heatmaps of enriched genes in denoted pathways based on KEGG from all groups [32][33][34] . Heatmaps were generated using NetworkAnalyst.

Discussion
Establishment of 3D culture models. We investigated the bilateral effect of the cellular architecture and milieu on the properties and phenotypic responses of the whole cellular microenvironment of tumor cells. Liver cancer microenvironment was modeled using 2D and 3D mono-and co-culture systems. The transcriptome of five models was profiled to determine the effects of each condition. Each one of the culturing conditions mimicked an aspect of the TME. Based on our clustering methods, group 5 clustered away all the other groups, specifically group 1, in both the PCA plot and hierarchical cluster heatmap (Fig. 2a-c). This is in line with the Venn diagram, where group 5 had the highest number of unique genes that are not shared with other groups (Fig. 2d). This clear distinction suggests a fundamental difference in gene expression between the proposed model of group 5 and all other culturing conditions, signifying a considerable change in gene expression as a result of 3D co-culture of cancer spheroids with fibroblasts.
Mimicking hypoxic conditions is not sufficient to mimic in vivo conditions. Treatment with hypoxia mimicking agent significantly upregulated genes associated with p53 and AP-1 networks (Table S1). AP-1 proteins (i.e. Jun, Fos, and ATF families) are activated by hypoxia and are frequently deregulated in cancer 50 . Deregulation of ATF-2 and its network have been implicated in liver development, regeneration, and cirrhosis 51 . In addition, ATF-2 has been reported to play a role in HCC resistance to sorafenib, in vivo 52 . The presented hypoxia mimicking conditions also induced several metabolism pathways including fructose, mannose, starch, and sucrose metabolism, as well as glycolysis and gluconeogenesis (Fig. 3). Nonetheless, without treatment with CoCl 2 , HepG2 spheroids (group 4 & 5) upregulated genes were associated with HIF1-α and HIF2-α ( Fig. 3a; Table S1). This is consistent with studies promoting cancer spheroids as candidate solid tumor model, as they recapitulate many aspects of in vivo tumors including hypoxia 53 . Hypoxia-mediated pathways in turn promote survival, angiogenesis, invasion, and metastasis 54 . Hypoxia also plays a role in lipid and steroid metabolism 55 . Steroids promote tumor immune evasion by suppressing T cell activation and subsequently effecting immune-therapy outcome [56][57][58] . This is consistent with induction of steroid synthesis and estrogen signaling pathways in group 4 and 5 (physiological hypoxia), but not group 2 and 3, suggesting that chemically mimicking hypoxia is not sufficient to recapitulate a more encompassing hypoxic condition ( Fig. 3; Table S1). Taken www.nature.com/scientificreports/ together, these results demonstrate how adjusting HepG2 cells from 2 to 3D culture introduces hypoxia and its associated hallmarks, which better represents in vivo cancer conditions.

Reconstitution of the dynamic bilateral interaction between cancer and stromal cells through
3D co-culture system. Studying the interaction between cancer cells, stromal fibroblasts, and their contribution to tumorigenesis remains a challenge. We, therefore, created a simplified system to mimic stroma-tumor interaction. HepG2 cells were either cultured with fibroblasts conditioned media only (group 3) to mimic oneway interaction or cultured with fibroblasts in a 3D-culture system (group 5) to mimic a two-way communication. When cultured with fibroblasts, HepG2 upregulated genes were associated with integrin cell surface interactions, in addition to urokinase-type plasminogen activator (uPA) and uPAR-mediated signaling. Integrin signaling pathway (ISP) regulates the interaction with the extracellular environment in response to intracellular cues 59 β1 integrins are overexpressed in many tumors, and blocking their signaling transduction reduces survival and tumorgenicity of many cancers, in 2D and 3D in vitro cultures, and in vivo [60][61][62][63] . ISP transduction has been reported to be promoted by other cell surface proteins, such as urokinase receptor (uPAR) 64 . Suppressing uPAR expression or disruption uPA/uPAR interaction have been reported to inhibit tumor progression and metastasis 65 . In addition, canonical pathways associated with the upregulated genes in group 3 and 5 included IL-6 mediated signaling. This is consistent with Integrin increased signaling, where enhanced IL-6/STAT3 signaling is promoted by β1-Integrin pathway 66,67 . IL-6 overexpression has been reported in many cancers, including HCC, where it is suggested to promote the transition of fibroblasts to CAFs 68 . Indeed, studies have shown CAFs as the main source of IL-6, promoting survival, migration, invasion, angiogenesis and stemness in colorectal, gastric, and liver cancer cells 69 . By co-culturing HepG2 cells with fibroblasts, we were able to recapitulate signaling pathways essential in tumor-stroma crosstalk, creating a more reliable model to study complex TME interactions. Recent models that incorporate fibroblasts in direct or in in-direct co-culture with HCC spheroids have reported a tumor promoting effects of fibroblasts or their conditioned media 70,71 .
Recapitulating signaling pathways in HCC through 3D co-culture system. Insulin/IGF signaling pathway is known to be activated in many cancers including HCC. Studies have shown its essential role in carcinogenesis and metastasis 72 . The insulin pathway was generally activated in our proposed 3D co-culture model (group 5) in comparison to the other culturing conditions. Many of the factors involved in insulin pathway-in addition to factors involved in hypoxia, angiogenesis, and MAPK signaling-were upregulated in our study (Figs. 3 and 4); which was confirmed on the transcriptome and secretome levels (Fig. 4). HK1 is dramatically upregulated in group 5 (Fig. 4a). HK1 shows involvement in glycolysis, HIF, and insulin signaling pathways (Fig. 3a,b), a role that has been reported in various studies; where HK1 contributes to glycolysis, proliferation, migration, and invasion of HCC 49,73,74 . Similarly, PFKFB3, a direct target of HIF1, is upregulated in all groups in comparison to group 1 and regulates glucose metabolism and promotes cancer progression and growth 75 .
Overexpression of PFKFB3 is associated with poor prognosis of HCC, and its inhibition resulted in suppression of HCC growth in vitro and in vivo 76 and reversed the in vitro sorafenib-resistance of HCC cells 77 . EGLN3, that mediates crosstalk between hypoxia and insulin signaling pathways 78,79 , was upregulated in HCC hypoxic settings 80,81 . Consistently, EGLN3 was a significant factor driving hypoxia and insulin pathways in our proposed model (group 5). Culturing HepG2 under 3D co-culture conditions also enhanced the expression of genes promoting angiogenesis, migration, and invasion such as SERPINE1, ETS1, MMP3, and MMP7. ETS1, MMP3, and MMP7 are only upregulated in group 5 in comparison to group 1. ETS1 is involved in upregulating hypoxiatarget genes such as MMP3, and MMP7 82,83 . Downregulation of ETS1 was reported to inhibit metastasis and invasion of liver cancer cell lines 84 . Similarly, Expression of MMP3 and MMP7 is correlated with enhanced metastatic phenotype, where their inhibition suppressed invasion and migration of HCC cells 18 . SERPINE1 was distinctively upregulated in all groups compared to group 1, yet, with the highest fold change exhibited in group 5. Increased expression of PAI-1 (encoded by SERPINE1) is correlated with aggressive cancers and poor prognosis, where it is also associated with migration, invasion, and angiogenesis in HCC tissue 85 . Culturing HepG2 under 3D co-culture conditions also downregulated the expression of genes involved in cell cycle regulation and survival, including E2F2, E2F7, and E2F8. E2F transcription factors were only differentially expressed in group 5, where they were found to be downregulated in comparison to group 1, as confirmed by qPCR (Fig. 4a). Downregulation of E2F2, E2F7, and E2F8 prevents cell cycle arrest and enhances clonogenic survival 86 . Taken together, these findings indicate that culturing HepG2 in a 3D co-culture system (group 5) successfully recapitulate many signaling pathways that are important in HCC in vivo.
Pathway complementation analysis highlights the involvement of regulatory miRNAs. The cross-talk between cancer cells and other cells in the TME is also partly mediated through expression of miRNA and/or release of miRNAs through extracellular vesicles (EVs) 87 . Gene-miRNA interaction network analysis of significant DEGs of group 5 revealed a number of predicted regulatory miRNAs (Fig. S5, full list in Table S2), many of which are signature of HCC deregulated miRNAs 48,73 . For instance, miR-335 had the highest number of connections in group 5 (Fig. S5a). Serum of HCC patients undergoing TACE was analyzed for circulating miRNAs level, where miR-335 level was associated with significantly poor prognosis 88 . Thus, EV miR-335 has been suggested as a novel therapeutic strategy as it was shown to be involved in proliferation and invasion both in vitro and in animal model 89 . Our data suggest an interesting role of EV-based communication that could be explored in the future using our novel 3D co-culture HCC model, which resembles a simplified setting of the TME. www.nature.com/scientificreports/ Secretome profiling in our 3D co-culture supports genes expression profiling and recapitulates in vivo signatures. Cellular communication between cancer cells and their surrounding is partially driven by secreted proteins and other soluble factors including cytokines, chemokines, and growth factors. Using antibody microarrays, we determined the levels of different factors and their binding proteins that play a role in insulin pathway, angiogenesis, and cytokine signaling (Fig. 4b-f). Our results show that secretion of insulin/ IGF pathway proteins IGF-2, IGFBP1, and IGFBP2 is increased by co-culturing HepG2 cells with fibroblasts under 3D co-culture conditions in comparison to mono-cultures (Fig. 4b,g). IGF-2 is upregulated in several tumors including HCC. Its overexpression was notably detected in HCC patient and was shown to induce liver tumor formation, proliferation and angiogenesis in mice 45 . IGF binding proteins (IGFBPs) are essential in the IGF signaling axis, where they bind with high affinity to IGF-1 and IGF-2, and have been reported in HCC patients 72,90 . Granulocyte-, macrophage-, and granulocyte-macrophage-colony-stimulating factors (G-, M-, GM-, respectively) secreted levels have all increased in 3D co-culture setting in comparison to mono-cultures (Fig. 4e,g). G-CSF, M-CSF, and GM-CSF have been shown to be involved in liver regeneration, fibrosis, angiogenesis, and initiation and progression of liver cancer 91 . Similarly, levels of factors involved in angiogenic pathway including growth factors, angiopoietins, and matrix metalloproteinases are increased by co-culturing HepG2 cells with fibroblasts under 3D culture conditions (Fig. 4c,f-h). Despite that the 3D model we are proposing herein is restricted to cancer cells and fibroblasts, present data interestingly shows that our 3D model (group 5 settings) increase the secretion of several cytokines known to orchestrate the cross-talk between the tumor and its immune TME (Fig. 4d,g), triggering pro-tumor inflammation and immunosuppression 92 . Our pathway complementation analysis also highlighted two of the enriched genes of group 5 as prognostic markers, namely IGF1R and EGLN3 (Fig. 4i,j). Both have been reported as indicators of poor HCC prognosis 81,93,94 .
In conclusion, we propose a novel 3D tissue culture model of liver cancer that better mimics in vivo settings, where cancer spheroids were co-cultured with 3D fibroblasts in a transwell system. Our model allows for studying the impact of co-culture on individual cell populations, and allows for studying paracrine methods of cellular communication such as extracellular vesicles. Our results showed that in vivo like conditions cannot be mimicked by simply growing cancer cells as spheroids, but by co-culturing them with 3D fibroblast with which they were able to crosstalk. This was evident by the upregulation of several pathways involved in HCC, and the dramatic jump in secreted factors and surface receptors by co-cultured cancer cells, many of which are also involved in tumor-stroma interactions. We have explored the aspects of HCC our proposed model mimic by combining transcriptome and small scale-secretome analysis, which could be expanded in the future to include proteome analysis. Different models have been suggested for studying different aspects of the TME, and have been described recently [95][96][97][98][99][100][101] . There are certain limitations to our suggested model such as absence of direct cell-cell contact as a way of cellular communication, which doesn't allow for visualizing the full effects of co-culture. However, our model allows for studying other indirect forms of communication and transfer of cellular cues such as secreted factors and extracellular vesicles. Due to the unavailability of liver derived immortalized fibroblast cell line, we opted for SV-80, a lung derived immortalized fibroblast cell line. Despite this being a limitation of our model, our results suggest that lung fibroblasts were still able to reproduce the enhanced phenotype after co-culture, suggesting the plasticity of the TME. Other limitations that can be addressed in future work are including more TME cellular components (i.e. immune cells) that would yield a more robust TME model, working with matched organ-derived cells, and inclusion of patient material to minimize artifacts that could be introduced by commercial cell lines [95][96][97][98][99][100][101] . Nonetheless, compared to the conventional 2D culture, the proposed model exhibits an increase in the expression of genes associated with development, progression, and poor prognosis of HCC. Our results correlated with a robust phenotype that better mirrors in vivo HCC, from gene expression to prognosis markers, and therefore, a more reliable platform for molecular understanding of HCC.

Data availability
The datasets generated during the current study are available from the corresponding author upon request. The resources and suppliers used in this study have been provided above.